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Abstract 



I— I Skin is one of the largest human organ systems whose primary purpose is the protection of deeper tissues. 

^^ As such, the skin must maintain a homeostatic balance in the face of many microenvironmental and 

f-H genetic perturbations. At its simplest, skin homeostasis is maintained by the balance between skin cell 

growth and death such that skin architecture is preserved. This study presents a hybrid multiscale 

mathematical model of normal skin (vSkin). The model focuses on key cellular and microenvironmental 

variables that regulate homeostatic interactions among keratinocytes, melanocytes and fibroblasts, key 

Q^ components of the skin. The model recapitulates normal skin structure, and is robust enough to withstand 

"""" physical as well as biochemical perturbations. Furthermore, the vSkin model revealed the important role 

I of the skin microenvironment in melanoma initiation and progression. Our experiments showed that 

^ dermal fibroblasts, which are an important source of growth factors in the skin, adopt a phenotype that 

"^j" facilitates cancer cell growth and invasion when they become senescent. Based on these experimental 

y]^ results, we incorporated senescent fibroblasts into vSkin model and showed that senescent fibroblasts 

transform the skin microenvironment and subsequently change the skin architecture by enhancing the 

growth and invasion of normal melanocytes as well as early stage melanoma cells. These predictions 

^T are consistent with our experimental results as well as clinical observations. Our co-culture experiments 

show that the senescent fibroblasts promote the growth and invasion of non-tumorigenic melanoma cells. 

We also observed increased proteolytic activity in stromal fields adjacent to melanoma lesions in human 

r* histology. This leads us to the conclusion that, senescent fibroblasts create a pro-oncogenic environment 

. J^ that synergizes with mutations to drive melanoma initiation and progression and should therefore be 

S^ considered as a potential future therapeutic target. This study suggests a potential link between aging 

H in the skin microenvironment and the development of melanocytic neoplasms. 

Author Summary 

Skin homeostasis depends upon the complex interplay of skin cells as well as interactions between cells and 
the microenvironment. Here, we generated a virtual skin (vSkin) model in order to test our hypothesis 
that dysregulation of cell-microenvironment interaction leads to aberrant skin structure and furthermore 
recapitulates pathologic conditions of the skin. To this end, we used the hybrid cellular automata method. 
The model couples a cellular automata that describes biological rules for cell types with partial differential 
equations that describe continuous microenvironmental factors. Our simulation recapitulated normal 
skin growth as well as self-repair. It also showed that microenvironmental factors produced by stromal 
cells (fibroblasts) play an important role in maintaining normal skin homeostasis. We find that the 



transdifferentiation of fibroblasts due to senescence changes the skin microenvironment to drive melanoma 
initiation and progression. Our model gives new insights into the processes of melanoma initiation and 
progression, and suggests a novel strategy for treatment. 

Introduction 

Skin is the largest organ of the body. It is a physical barrier that limits the flow of water and electrolytes, 
while providing protection from ultraviolet radiation, microorganisms and toxic substances. Human skin 
consists of three layers, the epidermis, dermis, and subcuits. The epidermis consists of keratinocytes, 
melanocytes, Langerhans cells and Merkel cells. The dermis is composed of connective tissue, fibroblasts, 
blood vessels and lymphatics. The subcutis is made up of loose connective tissue and insulating fat 
cells. Keratinocytes are the main constituent of the epidermis, composing around 95% of the total cell 
number 1 . Within the epidermis the keratinocytes form a stratified squamous epithelium in which the 
keratinocytes are tightly connected to each other through desmosomes. Like other epithelial tissue, the 
epidermis continuously renews itself. In the epidermis, cell proliferation occurs at the basal layer followed 
by a stepwise upward migration and maturation of the cells. As the keratinocytes mature they lose their 
nuclei and turn into keratin- filled corneocytes that form a 10-30 cell thick water resistant protective layer 
that gives the skin its critical barrier function. Protection of the keratinocytes from the DNA-damaging 
effects of ultraviolet radiation is mediated by the pigment producing cells, melanocytes, which locate to 
the basement membrane of the skin in a 1:5 ratio with the basal keratinocytes [sjls]. Upon exposure to 
UV radiation, melanocytes produce melanin which is transferred to the surrounding keratinocytes via an 
active transport process. Once taken up into keratinocytes, melanosomes orientate themselves over the 
nucleus in a protective "hat" that shields the nuclei of the basal keratinocytes from UV-induced DNA 
damage 4 . Keratinocytes in turn regulate the growth and behavior of melanocytes through a complex 
network of cell-cell adhesion proteins and secretory factors 3 . In addition to these interactions, the 
structure and strength of skin is also critically dependent upon the extracellular matrix (ECM) and the 
fibroblasts in the dermis. Dermal fibroblasts secrete ECM components and are responsible for mechanical 
strength of the dermis 5 as well as producing growth factors such as epidermal growth factors (EGF), 
transforming growth factor /3 (TGF^), and basic fibroblast growth factor (bFGF) 6 . These growth 
factors promote growth and facilitate the constant renewal of the epidermis. 

Melanoma is the most devastating form of skin cancer [tIIs] , arising from the malignant transforma- 
tion of melanocytes. The first phenotypic change in melanocytes is the disruption of growth controls and 
development of melanocytic nevi, a benign mole. At the molecular level, the growth is stimulated by the 
abnormal activation of the mitogen- activated protein kinase (MAPK) signaling pathway which can arise 
through activating mutations in the NRAS (15% of melanomas) or BRAF (50% of melanomas) (£-12 



Despite BRAF mutations being important for melanoma development, the majority of benign nevi also 
harbor BRAF mutations 13 and rarely undergo full malignant transformation. Instead, benign nevi 
cease proliferation and remain quiescent for decades, entering a state of oncogene-induced senescence 
(OIS) [MJflS]. This suggests that nevi must acquire additional intrinsic (molecular) or extrinsic (mi- 
croenvironmental) damages to free themselves from the OIS state and develop into full-blown malignant 
melanoma. There is some suggestion that activation of phosphoinositide 3-kinase (PI3K) /protein kinase 
B (PKB or AKT) due to loss of a tumor-suppressor gene, phosphatase and tensin homologue (PTEN) may 
be required for escape from senescence 16 . The suggested mechanism is that the inactivation of PTEN 
fails to attenuate the PI3K/AKT level, and consequently increases PI3K/AKT expression that facilitates 
cell proliferation and survival. However, the precise molecular events that underlie the transformation of 
nevus cells develop into melanoma are not well understood. 

Cancer is typically a disease of old-age, and there is increasing evidence that senescence within the 
stroma, particularly in the fibroblast compartment can drive tumor development. A number of reports 
have shown that senescent stromal fibroblasts stimulate premalignant and malignant epithelial cells to 



grow in cell culture and to form tumors in mice 17-20 . Mechanistically, this seems to involve the secretion 



of factors from senescent fibroblasts, such as Matrix metalloproteinase-3 (MMP-3) and Interleukin-6 (IL- 
6), that in turn remodel the microenvironment, alter epithelial differentiation, promote endothelial cell 
motility and stimulate cancer cell growth both in vitro and in vivo 2l][22 . It is well known that aged skin 



tends to harbor large numbers of senescent fibroblasts 23 , and there is a direct correlation between age 
and the expression of pl6^^^^^ (a cell cycle inhibitor) expression within the dermis 24 . In other words, 
concentrations of pi 6^^^^^ increase dramatically as tissue ages. Photo-damage following exposure to 
UV radiation can also increase the level of senescence within the skin through the effects of oxidative 
damage upon the telomeres shortening [25|[26] . Although the exact role of the dermal microenvironment 
in melanoma development has been little studied, there is evidence that the mildly hypoxic environment 
of the skin contributes to melanocyte transformation ^^. It is also known that fibroblasts contribute 
towards the survival of early-stage melanoma cells by secreting the growth factor IGF-1 and that stromal 
remodeling through TGF^ overexpression can enhance melanoma survival in mice 
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In order to gain a better insight into the mechanisms underlying normal skin homeostasis, and ul- 
timately the role of the stroma in melanoma initiation and progression, we developed a virtual skin 
model. Several computational studies of skin have recently been developed. Deterministic ordinary dif- 
ferential equations have been used to model immune signaling in human skin, providing a quantitative 
strategy that distinguishes healthy from pathologic infiammatory responses [29] . A similar method has 



been employed to model mouse embryonic melanoblast proliferation dynamics 30 . Another continuous 
modeling framework, a system of deterministic reaction-diffusion equations, has been applied to model 
tumor-immune interactions 31 . A continuum fiuid mechanics approach has been adopted to analyze the 



stability of the interface between epithelia and stroma 32 . A mixture theory has been used for structure 
stability analysis of melanoma growth [^ and microstructure patterning |34 . More discrete based mod- 
eling frameworks have also been applied to model skin. Agent based methods have been used to model 



epidermal keratinocyte behavior in normal skin 35-37 , epidermal homeostasis control 38 , keratinocyte 
originated skin disease [39], the interactions between keratinocytes and fibroblasts [40^, and melanocyte 



distribution in epidermis kll. However, none of these previous studies have considered normal skin 
homeostasis and its disruption as a route to melanoma development. 

In this study, we have developed a hybrid multiscale mathematical model that simulates normal 
skin homeostasis. We named the model virtual skin (vSkin). We employed a hybrid cellular automata 



(HCA) approach 42-45 . The general modeling approach is to derive a set of rules for discrete cells and 
couple this discrete cellular population with a suite of continuous microenvironment al factors. Previously, 
we have applied this method to investigate tumor-stroma interactions in prostate cancer progression 46] . 
In this study, we developed a minimal set of cell rules for melanocytes, keratinocytes, and fibroblasts, 
which determine each cell's life cycle. Then, we developed a cell interaction network that defines local 
interactions between cells and their microenvironment. These local interactions lead to the emergence 
of normal skin structure that recapitulates the in vivo structure of skin. Using our model, we observed 
that the vSkin model is robust enough to withstand perturbations. The vSkin model not only recovers 
from massive loss of its constituents but also finds an equilibrium state after super-physiological per- 
turbations of microenvironmental factors. The proposed model also reveals the possible involvement of 
these microenvironmental factors in conjunction with senescent stroma in driving melanoma initiation 
and progression. 

The paper is organized as follows. In the first section, we give an overview of the vSkin model 
development process. In the second section, we initially present results showing the stability of vSkin. 
Then, we discuss the robustness of the vSkin as well as application of the model to melanoma initiation 
via melanocyte mutation and fibroblast transdifferentiation. We conclude with a discussion focusing on 
the implications of the model predictions in terms of a novel strategy for treatment. 
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Figure 1. Human skin normal structure. A: H & E stained cross-section of human normal skin at 
20X showing the epidermis contains a basal layer, melanocytes and keratinocytes. The dermis contains 
fibroblasts and extracellular matrix. B: Model domain with its key cell types such as keratinocyte, 
melanocyte and fibroblasts. The density of extracellular matrix is considered and assumed to be 
continuous in the domain. The gray color represents the density of extracellular matrix at each node in 
the domain, i.e., darker gray represents denser matrix. The basement membrane separates epidermis 
and dermis layer in the model, and it is represented by continuously connected nodes where the density 
of matrix has reached its maximum (1.0). 

Material and methods 

Mathematical model 

We propose a systemic model of skin capable of growth control as well as self-repair. In this section, we 
describe the development of the vSkin model, which represents a cross-section of skin as shown in Figure 
[l] We first present key microenvironmental variables involved in normal skin structure and function. 
Then we move on to discuss the implementation of key cell phenotypes. 

Hybrid cellular automata model 

In order to keep the complexity of the model to a minimum, we only consider three key cell types critical 
for skin function, melanocytes, keratinocytes, and fibroblasts. These cells are defined as points on a two- 
dimensional grid that represents a cross-section of human skin (6 mm x 1.24 mm), as shown in Figure [lb. 
Each grid point may also be occupied by up to five microenvironmental variables, EGF, bFGF, TGF/3, 
MMPs, and basement membrane/ECM. The discrete cells are coupled with these microenvironmental 
concentrations to define the hybrid cellular automata (HCA) model. The cells on the lattice infiuence 
these microenvironmental concentration fields through production and consumption, but are in turn 
affected by the concentrations, as they serve as regulators of cell behavior. To better understand how these 
cells interact with one another and their microenvironmental variables we developed a cell interaction 
network (Figure |2]). This network shows the key interactions that are believed to be important in 
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Figure 2. Cell interaction network for normal skin. Green lines represent promoting events, 
while red lines imply inhibiting effects. Note, cell coloration is used in all simulations to represent the 
equivalent cell. In the model, a fibroblast is assumed to produce three growth factors. Epidermal 
growth factor (EGF), Basic fibroblast growth factor (bFGF) and Transforming growth factor (TGF^). 
EGF promotes growth of keratinocytes, whilst TGF/3 inhibits the growth. The growth of melanocyte is 
controlled by both bFGF and physical interaction (contact inhibition) with neighboring keratinocytes. 

maintaining normal skin structure and function. The network was derived through lengthy discussion 
with our experimental co-authors and an extensive literature review, detailed throughout the following 
sections. 



Microenvironmental factors 

Keratinocytes, melanocytes and fibroblasts require a large number of different growth factors and adhe- 
sion signals to promote and maintain their normal growth and maintenance. However, here we restrict 
ourselves to either one or two chemical variables for each cell type in the model for the sake of simplicity. 
For a keratinocyte, we consider EGF and TGF/3, because these two factors are known to be the most 
effective regulators for keratinocyte growth control 47-51 . We consider bFGF for the growth control 



of both melanocytes 52 and fibroblasts 53 . Since fibroblasts produce protease when they are acti- 
vated l5||6], we include MMP in the model as well. Finally, ECM concentration is incorporated, since 
ECM is required for structural support, and it separates the dermal and epidermal layers [lllSj in the 
form of basement membrane. The dynamics of the five microenvironmental variables are defined by a 
system of partial differential equations that describe how each of them evolves in space and time. Since 
these continuous variables interact with discrete cells, for clarity we present only the discretized version 
of each equation. 
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where a ceh type X can be keratinocyte (i^), melanocyte (M) or fibroblast (F). In Eq (fTl), [•] represents 
the concentration of a chemical at a node r] at time step t (e.g.,[G£;] = GE{t^r]x^r]y))^ [-]~^^ represents 
the concentration at the next time step, and St represents a time step. The operator A2 defines a 
two-dimensional discrete Laplacian, 

/v ./, X _ f{t, ^x + Sh, r]y) + f{t, rja: - Sh, r]y) + /(t, 77^, r]y + Sh) + /(t, r]^, r]y - Sh) - 4f{t, r^^^, rjy) 
A2f [t,r]^,r]y) = ^^— ^ , 

where Sh is the grid size. Equation (fTl) defines EGF as a diffusing chemical which is increased by a rate 
of aE or Pe when a grid node r] is occupied by a fibroblast (F) or a keratinocyte (K)^ respectively. The 
concentration is decreased if a node is occupied by a keratinocyte or melanocyte (M) at a rate of je and 
(e^ respectively, due to uptake of EGF. Lastly, the concentration has a natural decay rate of A per time 
step. Similarly, the dynamics of bFGF (G5) is modeled as follows. 
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where fibroblasts and keratinocytes produce bFGF at the rate of a^ and /35, respectively, and transformed 
fibroblasts and transformed melanocytes produce bFGF at the rate of a^ and /3g, respectively. All cell 
types bind to bFGF at a rate of 75, ^6, ^6, Q^ and ^^. Finally, TGF/3 (T/3) is produced by fibroblasts and 
keratinocytes at a rate of a^^ and /3t^ , respectively, diffuses at a rate Dt^^ , and binds to keratinocytes, 
melanocytes (or transformed melanocytes), fibroblasts (or abnormal fibroblasts) at rates of 7t^, Ct;3, ^T^g, 
respectively. 
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ECM (E) is produced by both fibroblasts or keratinocytes depending on the local concentration of 
ECM. In other words, keratinocytes and fibroblasts produce ECM with a rate of hZE and pe only if there 
is any loss compared to the initial density of epidermis {Eq) and dermis (Fi), respectively. However, 
transformed fibroblasts always produce ECM with a rate of p|;. Whenever fibroblasts are in contact with 
keratinocytes or melanocytes, fibroblasts can make much denser matrix which will create new basement 



membrane 54 . Basement membrane (BM) is defined as a continuously connected curve of grid points, 
each of whose ECM densities is maximal (taken to be 1 in non-dimensional units). ECM is degraded by 
MMP at a rate of /i. The governing equation for ECM is 
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where Eq represents initial density of epidermis and Ei is that of dermis. Note that when a fibroblast is 
in contact with either a melanocyte or a keratinocyte, we set Ei to be its maximal value 1.0 to model 
fibroblast's basement membrane generation. The function H is the Heaviside step function defined as 

Hi.)4\ ^'^ (5) 

II it X > 0. 

MMPs (P) are produced by abnormal fibroblasts, diffuse at a rate D^^ and are depleted as they degrade 
the ECM (E). Note that when melanocytes are transformed, they can produce MMPs as well. In the 
model, we assume that transformed melanocytes produce MMPs with a rate of /3^. The governing 
equation is 
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^^ DpA2[P]+ a^F^^ + P;m:^ - ^[P][E]. (6) 

All boundary conditions for all microenvironmental variables are no-fiux on the surface and the bottom 
of the domain. Periodic boundary conditions were imposed on the left and right sides of the domain. 

Cell phenotypes 

The behavior of each cell depends on its neighbors and the microenvironmental concentrations in which it 
resides. In the vSkin model, we consider the von Neumann neighborhood of range 1 (i.e., four orthogonal 
neighbors) for each grid point rj e ft, denoted by A/"^. In vSkin, each cell has three possible phenotypes: 
proliferative, migratory and dead. In this section, we describe how these phenotypes are implemented. 
Biological rules for each cell were abstracted from a literature review and are summarized as fiowcharts 
in Figure [3] For the sake of simplicity, each cell is allowed to execute only one phenotype, proliferative, 
migratory, or dead, at each cell time step {tc). In other words, if the cell chooses one phenotype at current 
time step, it exits current life cycle and waits until next time step. If a cell does not execute any of these 
three phenotypes, it remains as a quiescent cell. 

Cell- cell adhesion: In order to incorporate one of the important mechanical aspects of skin into the 
model, we define an adhesion level (0 < A^ < 1) at each node r] e ft explicitly by using cell-cell adhesive 
coefficients C^,^, where ^ is one of its orthogonal neighbors (i.e. ^ G 7V^). The adhesive coefficient between 
the same cell types is defined as 0.5. A value zero is assigned when there is no known adhesion mechanism 
between two cells, such as between a fibroblast and a keratinocyte. The coefficient for two melanocytes 
is based on current literature [3,55 . To model E-cadherin mediated adhesion between a keratinocyte 



and a melanocyte, we impose a higher value (1) between these cell types. To recapitulate the anchorage 
of a melanocyte to the basement membrane, we also assign a higher value (1). The adhesive coefficient 
C^^^ for all cases is summarized in Figure [4] A-B. By averaging the four adhesion coefficients C^^^ at each 
grid point, we define an adhesion level A^ for a cell at that grid point rj. 

Av=T. C.,,«/4.0. (7) 



i 



Produce 
EGF, bFGF & 
TGPP __ 



1 - 



Migrate 



^ Sufficient 
"^ bFGF 



/ 



Keratinocyte 

\ 

Produce 
bFGF&TGFp 

\ 

Migrate ^ 

\ 

Sufficient w 
EGF ^ 



H Proliferat^H 



\ 




Death 



Sufficient 
TGFp 



\ 



Sufficient 
bFGF 



Proliferate 




.Yes 
► No 



Figure 3. Flowcharts showing the cell life cycle for fibroblasts, keratinocytes and 
melanocytes. A: Fibroblasts produces bFGF, EGF and TGF^, and proliferate if there is enough 
bFGF. Low concentrations of bFGF trigger fibroblast death. B: Keratinocytes produces bFGF, and 
proliferate when there is enough EGF and a low concentration of TGF^. An insufficient level of EGF 
triggers keratinocyte death. C: Melanocyte start their life cycle by checking the melanin unit (the 
number of keratinocyte neighbors) and then decide to migrate or not. Sufficient bFGF stimulates 
melanocyte proliferation while insufficient bFGF triggers cell death. 



Migration: In vSkin, every cell has capacity to move to one of its orthogonal neighbors (Nj^) regulated 
by its cell-cell adhesion restriction (i.e., Pr(migration) ~ (1 — A^)), where A^ is a cell-cell adhesion level 
obtained from equation ^. Since the migration rules are slightly different for each cell type we will give 
a more detailed description in the following paragraphs. 

Both keratinocytes and fibroblasts move only if empty grid points exist in A/^. We first determine 
the cell-cell adhesion value (A^^). If there is only one empty neighbor in the neighborhood N^^ the cell 
moves to this point with probability (1 — A^). If there are more empty grid points, the new position for 
the cell is chosen randomly from these grid points. An example of keratinocyte migration is shown in 
Figure [4p-D. 

In addition to the cell-cell adhesion restriction, we use more complex rules for melanocyte migration 
based on the current literature l2l|3]. Melanocytes locate to the basement membrane in a 1:5 stable ratio 
with basal keratinocytes. Typically, one melanocyte is associated with approximately 36 keratinocyte 
neighbors. This symbolic structural relationship is known as the melanin unit. In the vSkin model, this 
melanin unit is modeled as a number of keratinocytes in the upper half of the Moore neighborhood of 
range 7 for each melanocyte. The value of 7 is set to 4 in vSkin. In vSkin, a melanocyte is moved to 
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Figure 4. Cell-cell adhesion and migration A: The cell-cell adhesion value A^ at each node is an 
average of its four orthogonal cell-cell adhesive coefficients {Cr^^^^ g, defined in B). C: An example of 
cell- adhesion level for a keratinocyte at a grid point rj. This keratinocyte happens to have one 
keratinocyte and one melanocyte neighbor, from ^ an adhesion value A^ is determined (A^ = 0.375). 
D: This keratinocyte can move to either the north or south grid point (empty) with the same 
probability of Pr = (1 — Ajj)/2 = 0.3125. E: Another example, where a melanocyte happens to have 
three neighbors (eastimelanocyte, south:keratinocyte, west: fibroblast) which gives a cell-adhesion value 
of 0.5. F: The probability of movement for this melanocyte is 0.5 (i.e. (1 — Arj)) and the cell can move 
to either its north or south grid point. However, since there is an ECM gradient ([£^] (south) > [E]{r]) 
but [E'] (north) < [£^](7^))the melanocyte is not allowed to move north and therefore can only move 
south. ECM gradient will always dictate the movement direction, provided one exists and their is 
sufficient space to move. 



one of its orthogonal neighbors with regards to its cell-cell adhesion level, provided that (i) the melanin 
unit of the melanocyte is disrupted, (ii) there is an ECM gradient towards one of its neighbors, and (iii) 
at least one of its neighbors is either empty or only occupied by a keratinocyte. In order to satisfy (i), 
a check for the number of keratinocyte neighbors is made every cell time step. Rule (ii) is implemented 



10 



because melanocytes are known to have a delicate cell-matrix interaction mechanism 56|[57 and migrate 



via haptotaxis. Based on rule (iii), a melanocyte is allowed to be moved to another grid point ^ in N^^ even 
if the point is occupied by a keratinocyte. This additional rule is implemented to recapitulate observed 



melanocyte high motile activity 56 . If there exist one neighbor (^) satisfying rule (ii) and (iii), the 
melanocyte is moved to the new position ^. If more than one such neighbors exist, the new position for 
the melanocyte is chosen randomly from the neighbors. An example of melanocyte migration is provided 
in Figure |4p-F. 

Proliferation: In the vSkin model, we assume that each cell has capacity for proliferation and will produce 
two daughter cells provided that (i) the cell specific growth factor is sufficient for cell growth and (ii) there 
is sufficient space surrounding the parent cell for the two daughter cells to occupy. For a keratinocyte 
to divide, the concentration of Ge has to be sufficient {Ge > Bj.) for its division. Then, a keratinocyte 
proliferates with a probability of Pr (keratinocyte division) ^ (1 — Tf^). A melanocyte proliferates if (i) 
the melanin unit is disrupted and (ii) the concentration of bFGF is sufficient for cell division (G5 > 
Bm)' A fibroblast is allowed to proliferate if the growth factor bFGF (G5) is greater than the threshold 
{Bf). In addition to the bFGF level, we also consider a space constraint that limits the number of 
fibroblasts in vSkin domain. A fibroblast can divide only if at most one fibroblast neighbor exists in its 
neighborhood. This is not an unreasonable assumption in that the main role of fibroblasts in the model 
is to supply microenvironmental factors and the number of fibroblasts and the production rate of the 
microenvironmental factors are scalable. In other words, we can first choose a fixed number of normal 
fibroblasts in the model and then scale the production rates (a^;, 0^5, ar^) accordingly. In order to satisfy 
(ii), we assumed that one daughter cell replaces the parent cell and the other daughter cell will move to 
any one of the parent cell's four empty orthogonal neighbors (A/"^). If more than one of the neighboring 
grid points is empty, then the new cell position is chosen randomly from these points. 

Death. For any cell to survive, it requires sufficient growth factors. We assume that the concentration has 
to drop to starvation levels, Sk for a keratinocyte and Sm for a melanocyte, before death can occur. For 
a normal fibroblast, we use a different constraint that incorporates both cell age and the concentration 
of bFGF. Since every fibroblast is continuously producing bFGF with the rate a^, it is never deprived of 
growth factor in normal skin conditions. We make another assumption that as a fibroblast ages, it will 
need more growth factors. In other words, the probability of fibroblast death is positively correlated with 
age but has a reciprocal relationship to the concentration of bFGF (Pr (fibroblast death) ~ [G5]/age). 
The space that dead cells occupy immediately becomes available to new cells. 

HCA Algorithm 

We now summarize the computational algorithm that integrates the discrete cells and the continuous 
microenvironmental variables. 

1. Define a rectangular grid {Q) that determines both microenvironmental concentrations and cell 
positions. 

2. Set initial concentration of EGF, bFGF, TGF/3 and ECM at each node r] eVt. 

3. Place fibroblasts, melanocytes and keratinocytes on the grid. 

4. Microenvironmental concentrations are solved using the governing equations (fTl) - ([6|. 

5. Cell-cell adhesion level A^ is determined for each rj G Q. 

6. The concentration of microenvironmental factors and adhesion level are coupled with each cell to 
determine its phenotype. 

7. Action associated with the determined phenotype is realized and the grid is updated accordingly. 
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8. Go to step 4. 

Parameterization 

This vSkin model inevitably contains a large number of parameters as it has three different cell types 
and five microenvironmental factors. Unfortunately, many of the model parameters were difficult, if not, 
impossible to obtain, and therefore need to be estimated. However, due to the interdependent nature 
of the variables, precise parameterization may not be as important, rather the interactions between the 
variables is the real driving force. Homeostasis is our primary goal and as such served as central driver 
for parameterization. We do not doubt that other parameter sets could achieve similar outcomes. Note 
that all parameters are non-dimensionalized f45l. 

The diffusion rate of EGF {De) is estimated using the relation given in 58 with the molecular 
weight 133.07 kDa |59J. The estimated value is De = 2.5 x 10~^ cm^/s. The production rate of 
EGF by a fibroblast has not been measured. Thus, we estimate it to be a 0.5 (dimensional value = 
1.2 fg day" cell" ) from multiple simulations varying this parameter within range of [0,1] to achieve 
normal growth (net growth ~ 0) of keratinocytes. Fibroblasts are the main growth factor supplier in 
skin 5,6 . Thus, we make an assumption that a keratinocyte produces EGF with a significantly smaller 
rate (0.005; non-dimensional). Growth factor consumption rates are difficult to measure, especially at 
the level of single cells. We consider a normal activity rate of growth factor receptors. For example, 
we use 0.01 for a normal activation rate of EGF receptors of keratinocytes, meaning that a keratinocyte 
is able to bind 1% of EGF at the grid point (per cell) per time step. The rate is set to be a smaller 
value (0.001) for melanocytes since EGF is known to stimulate keratinocyte growth mainly 47{[48] . 
Note that we varied these consumption rates within a range [0, 1] to achieve skin homeostasis. We 



estimate the diffusion rate of TGF/3 using the relation in 58 and the molecular weight 44.3112 kDa [59|. 



The estimated value is Dt^ = 5.8 x 10 ^ cm^/s. The production rate of TGF^ by a fibroblast is 



chosen to be 0.21333 fg day cell from in vitro study 60 . The rate for a keratinocyte is assumed 



to be significantly smaller than that of a fibroblast 61 . Since TGF^ is known to mainly regulate the 



growth of keratinocytes 49J-I51], we assumed that keratinocytes bind to TGF/3 the most. The rates for a 
melanocyte and a fibroblast were estimated to achieve normal skin homeostasis. The consumption rates 
of TGF/3 are estimated to be 0.2, 0.6, 0.1 for a melanocyte, a keratinocyte, and a fibroblast, respectively. 
Using the same method [58' with molecular weight 17.1531 kDa, the diffusion rate D^ is estimated to be 
1.2 X 10~^ cm^/s. The production rate of bFGF by a fibroblast is also estimated from multiple simulations 
varying these parameters within a range of [0,1] (non-dimensional). The estimated parameter is 0.05 (0.12 
fg day" cell" ). The rate for keratinocytes is assumed to be significantly smaller, since keratinocytes 
appear to produce significantly less bFGF than fibroblasts 62 . Since both melanocytes and fibroblasts 
use bFGF as a growth promoting signal and there is no evidence showing which cell binds to bFGF more, 
we use the same value (0.02) for the bFGF binding rate. The rate for keratinocytes is set to be significantly 
smaller (0.005). The diffusion rate of MMP is estimated to be Dp = 5.01 x 10"^ cm^s in 58jj59j. The 
production rate of MMP by fibroblasts is varied within a range 0.024 - 0.12 fg day" cell" . Estimates 
for the kinetic parameters A and /i are not available since they are very difficult to obtain experimentally. 
Thus, we set non-dimensional values that best produced homeostasis. We know that ECM density in 
the dermis is significantly higher than that in the epidermis (Figure fl]), we therefore set the epidermal 
ECM density to 0.2 and the dermal ECM density is 0.7. All of the parameter values described here are 
summarized in Table [H 

The cell cycle time is set it to be a typical value tc = 24 h. The size of the grid is set to be 300 x 62. 
The grid size Sh is set to be 6h = 20 /im (a typical cell diameter), and thus the dimensions of the skin 
slice we simulate are 6 mm x 1.24 mm. The time step for the microenvironmental variables is set to 
be St = 4.8 h. We assume that all cell types have the same volume 3.1416 x 10~^ cm^ with a radius 
10 /im, and therefore the base number of cells is set to be 3.1831 x 10^. The background density of 
ECM is set to be 5.0 x 10~^g cm^ |63 . The background concentration of microenvironmental variables is 
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Table 1. Model parameters 



Parameter 


Description 


Value Normalized value 


Dt, 


Diffusion rate of TGF^ 


5.8 X 10-« cm'Vs 58 59 0.0838 


OtTp 


Fibroblast TGF/3 production rate 


0.21333 fg day"^ceir^ 


60 


0.0873 


Pt, 


Keratinocyte TGF^ production rate 


0.0106 fg day cell"^ 0.004 


IT, 


Keratinocyte TGF^ binding rate 


0.6 


Ct, 


Melanocyte TGF^ binding rate 


0.2 


^T, 


Fibroblast TGF^ binding rate 


0.1 


De 


Diffusion rate of EGF 


2.5 X 10"*^ cmVs 58 59 0.0367 


Q-E 


Fibroblast EGF production rate 


1.2 fg day celP^ 0.5 


Pe 


Keratinocyte EGF production rate 


0.0012 fg day cell"^ 0.005 


IE 


Keratinocyte EGF binding rate 


0.01 


Ce 


Melanocyte EGF binding rate 


0.001 


Db 


Diffusion rate of bFGF 


1.2 X 10"^ cmVs 58 59 0.1708 


otb 


Fibroblast bFGF production rate 


0.12 fg day celP^ 0.05 


13b 


Keratinocyte bFGF production rate 


0.024 fg day cell"^ 0.01 


^b 


Fibroblast bFGF binding rate 


0.02 


lb 


Keratinocyte bFGF binding rate 


0.005 


a 


Melanocyte bFGF binding rate 


0.02 


A 


Decay rate of all growth factors 


0.01 


D, 


Diffusion rate of MMP 


5.01 X 10-8 cmVs |58, 59 0.0724 


KE 


Keratinocyte ECM production rate 


15.7 pg day-^ celH*"^ 0.1 


PE 


Fibroblast ECM production rate 


15.7 pg day-^ celF^ 0.1 


IJ' 


ECM decay rate 


0.001 


Eo 


Epidermis ECM density 


0.2 


El 


Dermis ECM density 


0.7 



assumed to be Go = 10 ng cm~^ 64 . We assume that the concentration of growth factors above which 
a ceh proUferates in normal skin is set to the non-dimensional threshold value 0.02 for Bm^ 0.02 for Bf 
and 0.01 for Bk- Since the concentrations below which cells die from starvation are also variables, we 
set a starvation threshold to be 1% (5'/c), 2% (5'm), and 0.1% of the background microenvironmental 
concentration Gq for keratinocyte, melanocyte, and abnormal fibroblast respectively. 

Experimental methods 

Cell culture 

Melanoma cells were kindly provided by Dr Meenhard Herlyn (The Wistar Institute, Philadelphia, PA) 
and were cultured in RPMI medium supplemented with 5% fetal calf serum (FCS). 



Co-culture growth 

Senescence was induced in primary fibroblasts following irradiation with 10 Gy followed by 7 days of 
recovery as described in Campisi and colleagues |21 . Senescent and pre-senescent fibroblasts were plated 
upon 10 cm tissue culture plates overnight before the addition of Adeno-GFP tagged melanoma cells (a 
gift from Dr Amer Beg, Moffitt Cancer Center). Images of the growing cells were taken after 24, 48 and 
72 hours with a Nikon-TSlOO inverted fiuorescence microscope. 
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Western blotting 

Cells were plated and grown overnight. To detect ADAM-9 levels, ADAM-9 antibody (ab) was used 
for immunoblotting (IB). For all IB experiments, proteins were denatured prior to separation on 6-18% 
Tris- Glycine gels. Proteins were transferred to polyvinylidene fluoride (PVDF) and blocked with 5% 
non-fat milk in Tris-Buffered Saline Tween-20 (TBST). Blots were incubated overnight in primary ab 
diluted according to the manufacturers instructions. Secondary ab's conjugated to the enzyme horseradish 
peroxidase (HRP) were used for detection by chemiluminescense. 

Zymography 

Cells were plated and allowed to grow overnight. To detect functional ADAM-9 activity, proteins were 
not denatured prior to separation upon a casein gel. Cleavage of casein was imaged using an Epson V300 
photo scanner. 

Three-dimensional spheroid assays 

WM793 melanoma spheroids were prepared using the liquid overlay method 65 . Spheroids were treated 
with serum free RPMI, pre-senescent fibroblast conditioned RPMI or senescent fibroblast conditioned 
RPMI. Pictures of the invading spheroids were taken using a Nikon-TSlOO inverted fiuorescence micro- 
scope. Images were converted to gray scale and the spheroid frames were outlined using Image J (Image 
Processing and Analysis in Java) by modifying the threshold. The spheroid frames were than inverted 
and absolute intensity of the invading cells was quantified using Adobe Photoshop. 

Human specimen procurement 

A representative tumor sample was obtained from a patient who underwent surgical resection for metastatic 
melanoma and was prospectively consented and accrued to an existing melanoma tissue procurement pro- 
tocol approved by the Moffitt Cancer Center Scientific Review Committee and The University of South 
Florida Institutional Review Board. After deparaffinization, the slides were stained for ADAM9 (Chemi- 
con) before a thorough washing and incubation with a Alexa-Fluor 647 secondary antibody. The slide 
was also stained with DAPI to indicate the nuclei and analyzed using a Leica confocal microscope. Areas 
of stroma, tumor and leading edge were delineated through examination of corresponding H & E stained 
slide by a dermatopathologist. 

Initialization 

Using the HCA algorithm and parameterization described above, we first ran an initial simulation to 
obtain the starting configuration of the domain (Figure [l] B) for all subsequent simulations. We follow 
the experimental steps in the in vitro 3D organotypic skin reconstruct experiment Figure [5] A-B, as if 
we develop a virtual organotypic skin reconstruct. This 3D organotypic culture system has served as a 
foundation for many basic science studies as well as a skin transplantation model, and it is known to 
be very stable and homeostatic (see review in \66 ). The initial simulation starts with a dermal layer 
mixed with fibroblasts that we simulate for two weeks until dermal fibroblasts produce enough growth 
factors and ECM for keratinocyte and melanocyte growth. Then a mixed population of keratinocytes 
and melanocytes are placed on the top of this matured dermis. As soon as keratinocytes are placed, 
they rapidly grow until they fill almost the entire domain. Soon after the most of keratinocytes become 
quiescent and start to die. Finally, the population becomes stabilized. The melanocyte population also 
rapidly increases initially but soon finds its equilibrium. This initial simulation shows the key interactions 
in the model are well balanced and manifest a normal skin homeostasis. Figure [5p shows how the skin 
structure develops over the period of a year as well as the concentrations of the microenvironemental 
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Figure 5. Organotypic skin and vSkin development. A: Scheme showing the generation of in 
vitro orgonotypic skin reconstructs. A stromal layer of collagen mixed with fibroblasts is placed on top 
of the acellular collagen, and the layer is incubated for 4 - 7 days to allow the fibroblasts to constrict 
the collagen. The epithelial cells are layered on top of the fibroblast-containing stromal layer, and then 
the reconstruct is air lifted to promote differentiation and stratification of the epithelial layer. B: 
Photomicrograph shows a mature organotypic culture of a normal skin. C: Snapshots of vSkin 
initialization. The steps of the simulation follow the same experimental steps used to develop 
organotypic skin reconstruct. vSkin development at time step 0, 14 and 365 days are shown in the first 
row, and the distributions of growth factors, age and force at the end point (t = 365 days) are shown in 
second and third rows. 



factors, the age of all cells and the adhesion field after one year. After a period of around 6-8 weeks, 
vSkin has already reached a homeostatic configuration. Reassuringly, the timescales of vSkin formation 
and the organotypic skin reconstruct are similar. We will utilize the resulting skin structure that emerges 
after a year as our initial condition in all subsequent simulations, since the skin structure that naturally 
emerged from this initial simulation is morphologically close to real skin. 



Skin fitness 

The degree of abnormality of vSkin was characterized by a defined metric, "skin fitness (/(t))." We 
quantify ratios of each cell compartment (ratio of keratinocyte to melanocyte and that of fibroblasts to 
melanocytes) at each time step {t) and compare with the ratios of the normal state (to). We also took 
changes of epidermal thickness into account to monitor skin fitness. Changes from time step to to time 
step t (to -^ t) in both compartmental ratios and epidermal thickness were weighted equally to determine 
the final skin fitness. The skin fitness is scored from -1 to +1, where +1 represents the maximal fitness 
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and a normal skin state, and is defined as follows, 



. Ar Ah 

/(t) = 1.0- (Wi—+W2 — 



(8) 



where r represents the sum of two ratios, the ratio of keratinocyte to melanocyte and that of fibroblast 
to melanocyte, A represents the change from time step to to t, h stands for a epidermal thickness, and 
wi and W2 are the weights. In this study we choose the same weight {i.e.^wi = W2 = 0.5). 
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Figure 6. Temporal evolution of each cell population. A: Mean populations dynamics of 
keratinocytes, melanocytes and fibroblasts (out of 50 realizations). Both keratinocyte and melanocyte 
populations initially increase after placement in the vSkin domain, but after two months the population 
is stabilizes. B: A sample of keratinocyte net growth over one year. Keratinocytes proliferate until the 
population reaches a maximum in the vSkin domain, then after about two months, keratinocytes start 
to commit cell death and the net growth begins to oscillate around zero. C: A sample of melanocyte net 
growth over one year. Initially, there are more births than deaths, but soon after the net growth 
oscillates around zero. D: A sample of fibroblast net growth over one year. 



Results 



The vSkin model maintains a stable homeostasis 

We already know that our vSkin model can recapitulate normal skin structure. As a typical CA method 
inevitably contains stochastic components, multiple realizations (> 50) were carried out for each simula- 
tion. Our simulations show that vSkin can reach a stable cell net growth rate and total cell number in the 
domain. Keratinocytes rapidly grow for the first two months and then reach a stable state (Figure [g] A 
and B), while both melanocytes and fibroblasts retain their initial population and the net growth remains 
stable (Figure [g] A, C, and D). These data indicate that our vSkin model rapidly achieves and maintains 
a stable skin homeostasis. 
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Figure 7. Wound healing. A: Positive correlation between healing time and wound size (depth) 
derived from a suite of wound healing simulations with a range to different wound depths (average of 
100 runs). The healing time was defined as the time it takes for the basement membrane to be 
continuously connected. The mean healing time increases as wound size increases. B: Three snapshots 
from a wound healing simulation, with a wound that was generated by removing a triangular shape slice 
with height of 0.96 mm and width of 1.86 mm. 



The vSkin model is robust to physical perturbation 

Once the stability of the vSkin model has been established, we next examined its robustness. Two 
different types of perturbations were applied to determine the robustness of our model system. First, we 
tested if our vSkin model can withstand massive loss of its constituents. To this end, a triangular shaped 
injury, which mimics an in vivo puncture wound to the skin, was created in the center of domain. As a 
new space was introduced into vSkin, cells nearby the space have an increased tendency to move toward 
the new space due to gradient in the force distribution. In the dermal layer, fibroblasts migrate into 
the injury site and start to produce growth factors (EGF, TGF/3, and bFGF) and extracellular matrix 
proteins. These newly produced growth factors promote growth of keratinocytes in epidermis, thereby 
starting to fill the void (healing procedure). A new basement membrane is generated by interactions 
between fibroblasts and either keratinocytes or melanocytes. The recovery time was defined as the time 
it took the wound site to form a complete basement membrane. 

We also changed the size of the wound by increasing injury depth, and quantified the relationship 
between healing time and injury size. The recovery time from wounding tends to increase with wound 
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size (see Figure [t] A). One hundred realizations were carried out and the mean and standard errors are re- 
ported. Interestingly, the variation among realizations increases with wound size. This is presumably due 
to the stochastic nature of the cell migration, thereby in basement membrane generation. In other words, 
the development of new basement membrane is dependent on chance interactions between fibroblasts and 
epidermal cells (keratinocytes or melanocytes). An example snapshot of the wound-healing process is 
depicted in Figure [7|B. These wound healing simulations indicate that our vSkin model is robust in the 
face of significant physical perturbation. Another implication is that the vSkin model can recapitulate 
some basic in vivo skin dynamics. 
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Figure 8. Growth factor perturbations and skin fitness. Simulation results from maximizing 
(1.0) the concentration of each growth factor over the entire domain during 3, 9, 15 and 21 days 
(averaged over 50 runs). After the perturbed skin reaches its new equilibrium, the skin fitness was 
measured using (|8|. The EGF and TGF^ impulses have a negligible effect on skin fitness. The bFGF 
and MMP have a greater impact on skin fitness, reducing it to 60% and 0%, respectively. However, 
whilst the skin structure is quite close to normal (see upper inset vSkin histology) after bFGF 
perturbation, after MMP over-expression, the skin has quite a different structure (see lower inset vSkin 
histology) characterized by a lower basement membrane, and the epidermal cells (keratinocytes and 
melanocytes) have populated the dermal layer. 



The vSkin model returns to equilibrium after biochemical perturbations 

The second perturbation involves manipulation of microenvironmental factors to determine if the vSkin 
system is robust enough to withstand super-physiological microenvironmental changes. Specifically, each 
growth factor was maximized in the domain during 3, 9, 15 and 21 days. The skin fitness function f{t) 



18 

in equation ([8| was used to quantify abnormality of vSkin after perturbation. Note that we measure skin 
fitness f{t) when vSkin reached a new equiUbrium. 

When EGF impulses were given to vSkin, keratinocytes responded very quickly. Keratinocytes pro- 
liferated until their population saturated the domain and then became quiescent. After most of EGF has 
either decayed or been consumed by keratinocytes and melanocytes, the quiescent keratinocytes deprived 
of EGF start to die. The resulting skin histology appears to be normal as it preserves normal epidermal 
thickness and the normal ratio of melanocytes to keratinocytes. In Figure [8J the blue line shows what 
little impact the EGF impulse made on vSkin, even for the longest perturbation time. 

The effects of TGF/3 impulses on vSkin were also negligible even though the growth of keratinocytes 
was inhibited by high TGF^ concentration. This is presumably due to fast turn over rate of keratinocytes. 
The light blue line in Figure IS] depicts how the skin fitness changes over different periods of TGF/3 
perturbation. 

Upon applying bFGF impulses, melanocytes proliferated faster and started to occupy more space than 
keratinocytes. The overgrowth of melanocytes resulted in inhibition of growth of keratinocytes due to a 
lack of space. When overgrown melanocytes have consumed the excess bFGF, they started to commit 
apoptosis. Unlike EGF or TGF/5 impulses, vSkin obtained a new equilibrium state with a slight increase 
of melanocyte number, particularly near the basement membrane. The overall fitness score is 0.6 as 
depicted by the orange line in Figure |8] 

Microenvironmental factors can transform normal skin to pathologic skin 

Unlike EGF, TGF/3, or bFGF, the overexpression of MMP had profound impact on vSkin fitness. Upon 
applying MMP, the extracellular matrix and the basement membrane begin to degrade. Destruction 
of the basement membrane increased the probability of melanocytes and keratinocytes to migrate into 
the dermal layer. In particular, melanocytes had a higher tendency to migrate into the dermal layer 
(a new microenvironment for them) where they can find more growth factors and loss of control from 
keratinocytes. This new microenvironment stimulates rapid growth of melanocytes. As a result of MMP 
overexpression, the vSkin became thinner with both keratinocytes and melanocytes moving into the 
dermal layer (downward). vSkin was able to recover from small changes induced by short term MMP 
overexpression (3 days), but fails to compensate for long-term MMP overexpression (9, 15, and 21 days), 
as shown in Figure [S] 

Collectively, the results indicate our vSkin model is able to return to its original homeostatic state 
or to find a new equilibrium in the face of significant microenvironmental deviations; however, this is 
dependent upon the perturbations only being to those variables present in normal skin (e.g., EGF, TGF^, 
and bFGF). When vSkin was disturbed by a factor not typically present in normal skin, its robustness 
is dependent on the duration of exposure. With long term exposure of MMP, vSkin transforms to a 
pathologic state. This implies that the regulation of microenvironmental factors contributes to vSkin 
transformation from a normal to an abnormal, pathologic state. 

Homeostatic disruption can lead to melanoma initiation 

Our in silico perturbations have shown that vSkin is a robust system which recapitulates normal skin 
form and function. These experiments have also shown that bFGF and MMP are important microen- 
vironmental regulators in skin homeostasis, particularly in melanocyte homeostasis. Melanocytes are 
the cell of origin for one of the most malignant forms of skin cancer, melanoma. A key to understand- 
ing melanoma initiation is determining how transformed melanocytes disrupt skin homeostasis. As we 
have readily shown in the previous section this homeostasis is crucially dependent on both keratinocyte 
regulation of melanocytes as well as microenvironmental regulation. Fibroblasts are the primary source 
for many of the microenvironmental factors that regulate skin homeostasis. Therefore, disruption or 
transformation of fibroblasts should also have a significant impact on skin homeostasis. Intriguingly, 
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Figure 9. Phenotype changes in fibroblasts. A: SA-/3-Galactosidase staining where the blue cehs 
are SA-/3-Galactosidase positive (indicating senescence). B: Automated inverted fluorescence images of 
pre-senescent (left) and senescent (right) fibroblasts. C: Micro-array analysis shows the fold increase in 
expression of multiple proteins in senescent fibroblasts compared to normal fibroblasts. Specific families 
of growth factors, proteases and extracellular matrix protein show increased expression. D: Senescent 
fibroblasts showed enhanced expression of ADAM9 proteins. E: Zymography shows senescent fibroblasts 
enhance matrix degrading capacity. 



the major characteristic of transformed melanocytes that distinguishes them from normal melanocytes is 
their ability to produce bFGF and MMP. Similarly we now know that senescent fibroblasts show altered 
patterns of ECM expression, growth factor release and protease activity 21 . 

To investigate this secretory phenotype further we undertook our own experimental investigation 
of the changes that occur in fibroblasts when they become senescent. Irradiation (10 Gy) of human 
primary skin fibroblasts led to their entry into a senescent-like state, as demonstrated by increased 
staining for /3-galactosidase (Figure [9]A.,B). The onset of senescence in the fibroblasts was associated with 
a marked change in their mRNA profiles with increases noted in ECM constituents, growth factors and 
proteases (Figure [9]C). Of these, ADAM9 is a protease thought to be critical for melanoma/fibroblast 
interactions 67 . Western blotting demonstrated that ADAM9 expression was highly upregulated in 
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Figure 10. Cell interaction network for abnormal skin. Key interactions among normal cells 
(keratinocyte, melanocyte, fibroblasts), abnormal cells (senescent fibroblast and transformed 
melanocyte) and microenvironmental variables represented as network. Note, cell coloration is used in 
all simulations to represent the equivalent cell. Green lines represent promoting events, while red lines 
imply inhibiting effects. A transformed melanocyte can produce its own growth factor (bFGF) and 
MMP, while senescent fibroblasts can produce bFGF, MMP and extracellular matrix proteins. 



senescent fibroblasts compared to non-senescent fibroblasts and that this was associated with increased 
matrix degrading capacity (Figure [9)l),E). Therefore fibroblasts produce protease, increased levels of 
bFGF, and extracellular matrix proteins as they become senescent. 

In order to explore the impact of modifying these cell types on homeostasis, we incorporated trans- 



formed melanocytes and senescent fibroblasts into the vSkin model. Figure 10 shows how these altered 
cell phenotypes modify the cell interaction network, now incorporating normal skin cells, transformed 
melanocytes and senescent fibroblasts. 



Senescent fibroblasts aid melanocyte proliferation and invasion 

To tease apart the relative contributions of the two abnormal cell types we initially look at each of 
them separately. First, we studied the role of senescent fibroblasts by examining how the degree of 
senescence in a subpopulation of the fibroblasts transforms normal skin structure (Figure 11). The 
degree of senescence specifically refers to how much bFGF and MMP the fibroblasts produce, so the most 
senescent phenotype produces maximal bFGF and MMP. For each parameter set (bFGF production rate, 
MMP production rate), we carried out multiple realizations and quantified skin fitness at each time step 
{f{t)) using the average number of each cell compartment and the average epidermal thickness from 
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Figure 11. Degree of senescence and its effects on melanocyte growth. Senescent fibroblasts 
were placed in the center of the dermal layer, among normal fibroblasts, only 10% were chosen to be 
senescent. We considered a range of senescence in terms of the production rate of bFGF and that of 
MMP. Specifically, the nondimensional parameter space for senescent fibroblasts was 
(al.ap) e {{x,y)\0.05 <x< 0.8, 0.0 < ?/ < 0.04}, where ag represents bFGF production rate by 
senescent fibroblasts and a^ indicates MMP production rate by senescent fibroblasts. Multiple 
realizations (50) with each parameter set (ag,ap were performed, and the mean skin fitness was 
quantified at each time step using Q. The first row shows the skin fitness landscapes after two and six 
years. A2 and A6 represent skin fitness obtained from simulations with mildly senescent fibroblasts 
(bFGF production rate is 0.05 and MMP production rate is 0.01). B2 and B6 show the fitness with 
highly senescent fibroblasts (bFGF production rate is 0.8 and MMP production rate is 0.040). Lower 
panels show representative vSkin histology for each of these points and highlights how different they are 



from normal skin (see figure 10 for a key of color coded cell types). 



these realizations. We observed that any degree of senescence for the fibroblasts significantly enhanced 
growth of melanocytes and promoted invasion into the dermis. From the skin fitness landscapes (upper 
panel. Figure 11) it is clear that the intensity of skin disruption is directly dependent on the degree of 



fibroblast senescence. When senescent fibroblasts have weak secretory phenotypes (i.e., close to normal), 
the fitness level of the skin they produce is almost one (i.e., normal). The skin structure, however, is 
not completely normal as there is some accumulation of matrix near the senescent fibroblasts (see virtual 
histology Figure 11, A2). When senescent fibroblasts have stronger secretory phenotypes {B2 and Bq) 
the melanocyte population expands rapidly and skin fitness deteriorates to -0.4 {B2). Interestingly, the 
fitness subsequently increases to 0.0 as time progresses (see the difference between B2 and Bq). This 
is due to the rapid growth of melanocytes initially (^2), driven by the senescent fibroblasts, but this 
excessive growth leads to deprivation of the growth factors for both fibroblasts and keratinocyes (see 
Figure 10 for the interdependence of cell types and growth factors) resulting in cell death. Whilst not 
cancer, excessive melanocyte production is a key feature of melanoma growth. It is important to note 
that the melanocytes here are normal and are only responding to the abnormal microenvironmental cues 
being produced by the senescent fibroblasts. However, if these melanocytes were to become transformed 
then this would be true melanoma, therefore, senescent fibroblasts may play a critical role in creating an 
environment ripe for melanoma initiation. 
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Figure 12. Degree of melanocyte malignancy and its effects on melanoma progression. The 

degree of melanocyte transformation was varied by increasing production rate of bFGF and MMP. The 
parameter range was given as (/3|,/3p) G {{x^y)\0.01 < x < 0.05, 0.0 <y< 0.04}, where /3p represents 
bFGF production rate and /3^ describes MMP production rate of a transformed melanocyte. The first 
row shows the mean skin fitness landscape after 0.25 and 0.5 years (averaged over 50 runs). Ao.25,0.5 
represents the fitness from simulations with minimally transformed melanocytes (bFGF and MMP 
production rates of 0.01), while Bo.25,0.5 indicates the fitness achieved with more aggressive melanoma 
cells (bFGF production was 0.05 and MMP was 0.04. Lower panels show representative vSkin histology 
for each of these points (see figure 10 for a key of color coded cell types). 



Transformation of melanocytes drives melanoma initiation 

So far, our vSkin simulations have shown that senescent fibroblasts are ideal for disrupting the normal 
skin environment to enhance melanocyte proliferation and invasion. However, we know that melanocytes 
are significantly transformed in melanoma and so we also need to consider how such mutations will affect 
normal skin structure. We therefore integrated phenotypes of mutant melanocytes into the vSkin model 
in order to study the role of melanocyte transformation as a driver of melanoma initiation. Note that we 
assume all fibroblasts and keratinocytes are normal for these simulations. 

Just as we considered a spectrum of senescence in the fibroblast population, we examine a range of 
melanocyte transformation. Minimally transformed melanocytes have only lost growth inhibition and 
will proliferate even in the absence of bFGF. As malignancy increases, melanocytes gain the ability to 
produce bFGF and MMP at ever increasing levels. Therefore degree of malignancy specifically refers 
to how much bFGF and MMP the melanocytes produce, so the most malignant phenotype produces 
maximal bFGF and MMP. Interestingly, minimally mutated melanocytes (i.e. those that have only lost 
proliferative control) remained largely quiescent due to competition for space and growth factors as well 
as their inability to invade (results not shown). However, melanocytes with a greater mutational load 
proliferated rapidly and invaded taking over the dermal layer of the skin, as shown in Figure [12| and 



quantified by the change in fitness landscape at 0.25 years and 0.5 years. Since the more transformed 
melanocytes have autocrine signaling, they don't suffer from growth factor deficiency even with loss of 
fibroblasts. With any degree of malignancy, mutant melanocytes significantly lower skin fitness and do so 
in a very short time frame, contrast the 0.5 year landscape of Figure [12] to the 6 year landscape of FigurepT] 
The effects become more dominant as the degree of transformation increases (Ao.25,0.5 -^ ^0.25,0.5)- The 
mutant melanocytes disintegrate skin structure and completely take over as shown in a sample snapshot 
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at 5o.5- This vSkin transformation recapitulates a type of melanoma which progresses very rapidly. 




Virtual histology 




Figure 13. The synergistic effects of senescent fibroblasts and minimally transformed 
melanocytes in driving melanoma development. The degree of senescence was varied as in the 
Figure 11 The transformed melanocytes (blue) are assumed to be only minimal transformed, i.e., loss 



of growth inhibition. Mildly senescent fibroblasts A2,6 still manage to produce a small melanoma nodule 
(lower panel) but only impact skin fitness a little (upper panel). However, increasing the fibroblast 
senescence significantly impacts melanoma progression, destroying the skin (virtual histology, B2,6) and 
reducing the fitness all the way to zero (upper panel). See figure 10 for a key of color coded cell types. 



Transformed melanocytes exploit senescent fibroblasts to drive melanoma progression 

The next logical step was to integrate both transformed melanocytes and senescent fibroblasts together 
in the vSkin model, to investigate how interactions between these abnormal populations drive melanoma 
initiation and progression. In order to highlight the impact of these interactions we only considered 
minimally transformed melanocytes in the following simulation. It is worth noting that these minimally 
transformed melanocytes have only lost control over proliferation and are incapable of producing cancer 
in the vSkin model independently since subsequent transformation would be needed for that to occur 
(i.e. resulting in increased bFGF and MMP production). 

Figure [13] highlights how senescent fibroblasts help transformed melanocytes proliferate more rapidly 
and invade into the dermis. Even when minimally transformed melanocytes are combined with the 
least senescent fibroblasts in vSkin we see a positive impact on the melanoma growth, with a nodule of 
melanocytes (Figure flsK ) forming after 6 years. Increasing the degree of senescence in the fibroblast 
population results in a more striking effect with the transformed melanocytes taking over the dermal 
layer, especially those nearby senescent fibroblasts (Figure flSb). From these results we can conclude that 
strong secretory phenotypes (e.g. the fibroblasts in Figure |13| ^2,6) can transform essentially normal 
skin (with the exception of the initial mutant melanocyte) to a pathologic state (see density graph 52,6)- 
Therefore, senescent fibroblasts have the potential to propagate and maintain melanoma development, 
provided a minimally transformed melanocyte population already exists. 



24 



A WM 35 cell growth 

pre-senescent senescent 

24h| 



48h| 





* • 











72h 



B WM 793 cell invasion (48 h) 
control 





degree of invasion 
(A.U.) 



I 



pre-senescent senescent 



Figure 14. Senescent fibroblasts enhance melanoma cell growth and invasion. A: Senescent 
fibroblasts aid ttie growtti of early stage (non-tumorigenic) melanoma cells. Equal numbers of 
GFP-tagged (green) WM35 melanoma cells were seeded on top of either primary skin human fibroblasts 
(pre-senescent) or fibroblasts treated with radiation and allowed to undergo senescence (senescent) and 
allowed to grow for 72 hrs. Graph shows the mean (of three independent experiments) relative growth 
rates (relative to WM5 cells grown alone) where melanoma cells were plated on top of pre-senescent or 
senescent fibroblasts. B: Poorly invasive melanoma cells (WM793) were grown as spheroids alone, in 
co-culture with pre-senescent fibroblasts or in co-culture with senescent fibroblasts. The spheroids were 
then implanted in collagen and allowed to invade. Mean increases in invasion were calculated using 
Image J. 



In vitro experiments show senescent fibroblasts enhance melanocyte/melanoma 
growth and invasion 

So far, we have presented vSkin model predictions that senescent fibroblasts change the skin microenvi- 
ronment and produce an aberrant skin structure that results from the overgrowth and invasion of normal 
melanocytes and early stage melanoma cells. In this section, we compare the model predictions with an 
in vitro co-culture assay. 



25 



We explored the prediction that senescent fibroblasts enhanced the growth and invasion of minimally 
transformed melanocytes in a series of in vitro fibroblast co-culture models. It was noted that plating early 
stage WM35 melanoma cells (which are non-tumorigenic in mice) on top of senescent fibroblasts enhanced 



their growth compared to non-senescent fibroblasts (Figure 14 A,B). To study invasion in a tissue-like 
context, we next generated spheroids in which poorly invasive melanoma cells (WM793) were co-cultured 
with either non-senescent or senescent human skin fibroblasts in a 3D collagen gel (Figure Il4]C,D). It was 
noted that the WM793 cells invaded the collagen minimally when grown alone and that their invasive 
behavior was markedly increased when co-cultured with the senescent fibroblasts. Although not direct 
comparison with the vSkin results, our co-culture cell growth assay highlights one of predictions that the 
vSkin model makes, which is that senescent fibroblasts promote the growth of early stage melanoma cells. 
The potential clinical relevance of our model predictions and experimental findings was suggested by the 
observation that human melanoma specimens stained positively for the protease ADAM9 at the leading 



edge where the melanoma cells and fibroblasts interact (Figure 15). In contrast, little ADAM9 staining 



was noted in the tumor core where stromal infiltration was lacking. 

Discussion 

The regulation of skin homeostasis involves a complex interplay between different cell types and their 
microenvironment. To better understand how these interactions produce skin homeostasis, we developed 
a mathematical model that incorporates the key cellular (keratinocytes, melanocytes, and fibroblasts) 
and microenvironmental (TGF^, EGF, bFGF and ECM) components of normal skin. Our virtual skin 
model (vSkin) successfully recapitulated normal skin structure, skin homeostasis and wound healing 
responses. Surprisingly the homeostasis that emerged from this model was remarkably robust, being 
able to recover from both physical and biochemical perturbations. However, the recovery was altered by 
both the perturbation type and period. With short-term perturbations inducing increased population 
growth, followed by a return to stable homeostasis. More persistent microenvironmental perturbations 
(specifically with bFGF) led to an altered but stable homeostatic state characterized by an increased level 
of melanocytes. One key perturbation that the system wasn't robust to was over-production of MMP 
and even then it was only when MMP levels were elevated for significant periods of time. Although, 
this isn't too surprising as MMPs are not normally produced in normal skin. These results indicate the 
importance of microenvironmental factors in maintaining skin homeostasis and highlight both bFGF and 
MMP as key factors in disrupting this homeostasis. 

The dialogue between our three cell types, keratinocytes, melanocytes, and fibroblasts is primarily 
mediated by growth factors produced by fibroblasts and cell-cell contact. One of the primary reasons for 
incorporating fibroblasts into vSkin was to understand their role in normal skin homeostasis, disruption 
and cancer initiation. An important experimental observation that motivated this line of inquiry was 
that when fibroblasts take on a secretory phenotype (such as that seen when undergoing senescence) they 
produce increased levels of bFGF and matrix degrading proteases (Figure [9]), the very factors that normal 
skin seems to be most sensitive to. By introducing a population of senescent fibroblasts into vSkin we 
made the unexpected prediction that such secretory cells may play a central role in the initiation and 
progression of melanoma. The least senescent fibroblasts caused a temporary disruption of homeostasis 
but as we increased the level of senescence (i.e. increasing the levels of bFGF and MMP production) we 
were able to permanently disrupt the skin structure and produced melanocytic hyperplasia i.e. a mole 
like structure. However, if the senescent fibroblasts were to be removed, normal skin homeostasis would 
be reintroduced since the melanocytes are normal and only reacting to the altered microenvironment 
created by the secretory cells. 

Our vSkin model is the first model showing all the steps in skin carcinogenesis, from normal skin fol- 
lowed by melanocytic hyperplasia, and ending with melanoma. The activation of stroma due to senescence 
induces excessive growth of normal melanocytes, resulting in the development of melanocytic hyperpla- 
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tumor edge 



tumor core 




Figure 15. Melanoma tissue staining. Representative H&E- stained section of a nodule of 
metastatic melanoma (upper panel), stained with anti-ADAM9 (lower panel). The yellow line 
demarcates the tumor- stroma interface. There is high ADAM9 expression at the tumor /host interface 
in the melanoma specimen but weaker, more diffuse expression in the tumor core. 
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A Melanoma progression 




B Anti-MMP therapy 



— minimally mutated melanocytes 
— weak (+bFGF+MMP) senescent fibroblasts 

— moderate (++bFGF+MMP) senescent fibroblasts 
— extreme (++bFGF++MMP) senescent fibroblasts 

— bFGF producing melanoma cells 

— bFGF & MMP producing melanoma cells 
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Figure 16. Multiple routes of Melanoma progression and a potential Anti-MMP therapy. 

The total population of transformed melanocytes was quantified for all untreated (upper right) and 
Anti-MMP treated (bottom right) cases as a measure of tumor burden. Representative vSkin histology 
for each of these cases is also shown (left panels), see figure [lO] for a key of color coded cell types. A: 
The population of minimally mutated melanocytes (barely visible green line) remains constant. The 
most mutated melanocytes produce both bFGF and MMP and progress the most rapidly ((a), dark blue 
line). The melanocytes which produce bFGF alone, form melanoma in situ where the melanoma cells 
are constrained to the epidermis ((b), light blue line). The most senescent fibroblasts ((c), dark red line) 
help minimally mutated melanocytes to proliferate rapidly and progress to full melanoma. Senescent 
fibroblasts with weaker secretory phenotypes lead to a slower melanocyte proliferation and invaaion 
((d), red line and (e), orange line). B: Melanoma progression after anti-MMP treatment. Anti-MMP 
therapy effects the growth of all populations but is more effective for melanoma progression driven by 
senescent fibroblasts ((h)-(j)) than melanocyte transformation driven melanomas ((f) and (g)). 
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sia (Figure 11). The vSkin model also recapitulated melanoma development driven only by melanocyte 



transformation (Figure 12). In addition, the vSkin model proposes alternative routes to melanoma ini- 



tiation, by showing that melanoma can occur as a result of co-operation between mutant melanocytes 



and senescent fibroblasts (Figure 13). Given that most known driver mutations in melanoma are associ- 
ated with entry into oncogene-induced senescence (OIS), our model suggests the possibility that signals 
from the senescent stroma may aid melanocyte transformation by allowing them to exit OIS. There is 
already good evidence that increased PI3K/AKT signaling can overcome BRAF V600E-induced OIS in 
melanocytes 16 and significantly most of the growth factors and altered ECM interactions we describe 
here are known to exert their effects through the PI3K/AKT pathway 68 . 

Biochemical analysis shows that the alterations in senescent fibroblasts included increased expres- 
sion of multiple growth factors and proteases. These alterations then stimulated the growth of normal 
melanocytes and early stage melanoma cells. The increased level of protease in senescent fibroblasts 
enhanced the invasiveness of late stage melanoma cells. Among many proteases, ADAM9 was found to 
be a putative mediator of matrix remodeling in senescent dermal fibroblasts. Importantly, the analysis of 
ADAM9 expression in melanoma samples from patients revealed that the expression was increased at the 



host-tumor interface, but not within the tumor core (Figure 15). Our findings mirror previous studies 
implicating a role for ADAM9 expression in melanoma cell/fibroblast interactions at the invasive front 
that is required for MMP-1 and MMP-2 mediated matrix degradation 67 . Whilst not conclusive, the 
ADAM9 staining implicates the stroma is producing matrix degrading enzymes. This clinical observation 
along with all of the experimental results we presented were consistent with our model predictions that 
senescent fibroblasts contribute to melanoma initiation and early melanoma progression. 

Our vSkin model proposes multiple routes to melanoma initiation and progression, summarized in the 



upper panel of Figure 16 The obvious and rapid route shown by the blue lines is due to the accumulation 



of mutations in melanocytes to drive melanoma development. The alternative route involves cooperation 



of mutant melanocytes with senescent fibroblasts. Figure 16 (red lines). When senescent fibroblasts 
stimulate and maintain the excessive growth of minimally transformed melanocytes we observe the slowest 
melanoma growth. As senescent fibroblasts adopt stronger secretory phenotypes, melanoma progresses 
more rapidly. 

If our predictions stand up to future validation, the fibroblasts could represent an important target for 
melanoma chemoprevention. The vSkin transformation to melanoma suggests a possible novel treatment, 
by negating the infiuence malignant phenotypes it may be possible to normalize the skin microenvironment 
and return (or maintain) skin homeostasis. Since proteases are one of key contributing factors driving 
melanoma initiation and progression, both from the senescent fibroblast and transformed melanocyte 
perspective, anti-MMP therapy might be a good treatment option. Indeed, there is already evidence 
that mutations in MMP-8, leading to an inactivation of protease activity that shifts the pattern of 
matrix degradation and growth factor availability, is an important event in the genesis of ~23% of all 



melanomas 69 . As a purely hypothetical exercise we ran a suite of simulations to test this idea, almost 
as if we applied anti-MMP cream to the vSkin model such that MMP expression was blocked completely. 
Specifically, the MMP level is completely knocked down to zero uniformly in the whole domain, when the 
size of initial melanocyte population had doubled. This knock down is maintained until the simulation is 
finished. Our in silico anti-MMP treatment experiments show that the treatment was far more effective 
for the stroma dependent melanoma in that the population growth rates were significantly decreased (see 
Figure [16] (lower panel)). 

This study emphasizes the importance of normal skin homeostasis and its disruption in melanoma ini- 
tiation and progression. The possible role of MMP inhibition as a chemoprevention strategy for melanoma 
will require extensive future experimental validation. MMPs and other related proteases constitute a large 
family of enzymes with distinct pro and anti-oncogenic activities. The effective use of MMP blocking 
drugs is likely to require an exquisite sensitivity profile that current inhibitors do not yet offer. Our in- 
tegrated approach highlighted two key unifying factors in disrupting skin homeostasis, MMP and bFGF 
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production. Whether these factors were produced by transformed melanocytes or senescent fibroblasts led 
to similar skin degeneration. The most important implication is that we cannot only consider melanoma 
as a transformation of the melanocyte population but must also consider the fibroblast population. It is 
perhaps no accident then that melanoma really is a disease of the aged [7Q||72] , where fibroblast senes- 
cence may naturally occur as we age and aid minor melanocyte mutations (due to a lifetime of skin sun 
exposure) to become full blown melanoma. 
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